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Recent development of theoretical method for dissipative quantum transport have achieved no¬ 
table progresses in the weak or strong electron-phonon coupling regime. However, a generalized 
theory for dissipative quantum transport at arbitrary parameter regime is not figured out until now. 
In this work, a variational method for dissipative quantum transport at arbitrary electron-phonon 
coupling regime is developed by employing variational polaron theory. The optimal polaron trans¬ 
formation is determined by the optimization of the Feynman-Bogoliubov upper bound of free energy. 
The free energy minimization ends up with an optimal mean-field Hamiltonian and a minimal inter¬ 
action Hamiltonian. Hence, second-order perturbation can be applied to the transformed system, 
resulting an accurate and efficient method for the treatment of dissipative quantum transport with 
arbitrary electron-phonon coupling strength. Numerical benchmark calculation on an single site 
model with coupling to one phonon mode is also presented. 


PACS numbers: 73.23.-b,71.38.-k,72.10.Bg,73.63.-b 

I. INTRODUCTION 

Quantum transport and energy dissipation in nanos¬ 
tructures have been of great interest, which requires the 
study and understanding of nonequilibrium phenomena 
of open quantum systems. 1-6 In many previous stud¬ 
ies on dissipative quantum transport, the coupling be¬ 
tween the system and the phonon bath can be small. In 
this case, second-order perturbation (2PT) can be em¬ 
ployed, leading to a set of equation of motions (EOM). 6 
However, there are many open quantum systems with 
strong electron-phonon coupling. For instance, the cou¬ 
pling between system and bath in the energy transfer 
process in photosynthetic complexes is comparable to 
the electronic coupling; 7-9 Strong electron-phonon cou¬ 
pling is also observed in resonant tunneling in molecular 
transport junctions. 10,11 Even though there are several 
non-perturbative methods to evaluate the dynamics nu¬ 
merically exactly, such as hierarchical equations of mo¬ 
tion (HEOM), 12,13 quasi-adiabatic propagator path inte¬ 
gral (QUAPI) 14 and multiconfiguration time-dependent 
Hartree (MCTDH) approach, 15-17 these methods are 
computational demanding and not trivial to be imple¬ 
mented. 

There is another group of methods to deal with strong 
electron-phonon coupling, which makes use of polaron 
transformation. A polaron transformed second-order 
master equation for open quantum system has been de¬ 
veloped accordingly. 18-21 The polaron theory has also 
been employed to study the quantum transport in the 
strong coupling regime, 22-28 in which strong phonon 
blockade effect is observed. The polaron transforma¬ 
tion extends the validity of perturbative treatment of 
electron-phonon interaction to strong coupling regime, 
provided that the electronic couplings are small com¬ 
pared to the typical electron-phonon coupling and bath 
frequency. However, when this condition is not fulfilled, 
the polaron theory performs worse than the standard per¬ 


turbative approaches as the polaron cannot follow the 
system motion accurately. 

In order to overcome this shortcoming of polaron trans¬ 
formation based approaches, the variational method has 
been developed as a generalized polaron transforma¬ 
tion. 29-33 In contrast with full polaron transformation, 
the variational polaron approach searches for an opti¬ 
mal amount of transformation by employing variational 
principle, which is implemented by minimization of the 
Feynman-Bogolubov upper bound of free energy. 32 The 
optimal transformation is determined by the properties 
of the bath and system-bath coupling. Thanks to the 
variational principle, the variational polaron method is 
able to interpolate between the weak and strong coupling 
regimes and to capture the dynamic behaviour at arbi¬ 
trary parameter space. Moreover, compared to the many- 
body approaches, variational polaron method is compu¬ 
tationally economic. 

Inspired by the recently developed variational polaron 
master equation, 29,30 variational polaron method for dis¬ 
sipative quantum transport is developed in this work, 
which aims at handling the dissipative quantum trans¬ 
port at arbitrary electron-phonon coupling strength. The 
manuscript is organized as follows. Sec. II introduces the 
variational polaron theory for quantum transport and 
minimization method for determining the optimal pa¬ 
rameters for variational transformation. Starting from 
the variational polaron transformation, the transformed 
Hamiltonian is divided into two parts: a mean-field part 
and interacting part. Free energy optimization ends up 
with an optimal mean-field Hamiltonian and a small in¬ 
teracting perturbation. Then 2PT is applied to establish 
the quantum transport theory. The variational polaron 
theory for quantum transport is applied to study the 
quantum transport through a single site model with dif¬ 
ferent electron-phonon coupling strength, which is given 
in Sec. IV. Finally, Sec.V summarizes the method devel¬ 
oped in this work. 
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II. VARIATIONAL POLARON THEORY FOR 

QUANTUM TRANSPORT 

A. Model Hamiltonian 

Considering a system sandwiched between two leads. 
The electrons can transfer from one lead to the other 
through the system driven by bias voltage. The electrons 
can be scattered inelastically by phonons when trans¬ 
port through the system. The corresponding Hamilto¬ 
nian reads 

H = H s + Y / [ H a + H s a } + H B + H S B, (1) 

a 

where H$ is the Hamiltonian of the system; Hb is the 
Hamiltonian of phonon bath; H a is the Hamiltonian of 
lead a. The lead serves as electronic bath; H$b an d 
Hga are the Hamiltonians of system-phonon coupling 
and system-lead interaction, respectively. Expressions of 
those Hamiltonians are 

Hs =^e„4c„, 

n 

Hql — ^ ^ ^k a ^kr* 

k a 

Hsa. — ^ ^ ”1” H.C.], 

n,k a 

Hb ^ ^n,kbl l j-.bn,k 

n,k 

Hsb = ^ ' ffn,fc c n c w(^ra.fc 4" (2) 

n,k 

Where e n denotes the single-particle energy of the elec¬ 
tronic state n in the device and cI (c ra ) is the correspond¬ 
ing creation (annihilation) operator. Similarly, 6k a and 
4 (c ka ) represent the single-particle energy and creation 
(annihilation) operator of fcth electronic state in lead a, 
respectively. Vk aU is the system-lead coupling strength. 
The phonon creation and annihilation operators are rep¬ 
resented by 4 k and b n ,/.. And uj nt k and g n ^ are the 
corresponding the phonon energy and electron-phonon 
coupling constant. 

In the absence of electron-phonon coupling, an ex¬ 
act time-dependent quantum transport theory has been 
developed recently. 34-38 In presence of electron-phonon 
coupling, electron has the probability of being scattered 
inelastically by phonon. Previous theoretical studies of 
dissipative quantum transport usually focus on two dif¬ 
ferent regimes: weak or strong electron-phonon coupling. 
In the weak coupling regime, 2PT (or lowest order expan¬ 
sion) is widely employed. While in the strong electron- 
phonon coupling limit, the electron is localized by phonon 
scattering. At this regime, polaron theory is suitable for 
describing the phenomena. In this manuscript, a gen¬ 
eralized theory for dissipative quantum transport with 
arbitrary electron-phonon coupling is developed by ex¬ 
tending the polaron theory to its variational version. 


B. Variational polaron transformation 


Similar to the polaron theory, the variational polaron 
theory starts from an unitary transformation generated 
by the operator 


U = exp 


-^n,fc4 c «(4,fc 

n,k 


bn,k) ? 


( 3 ) 


where A Ut k = fn,k/^n,k- The unitary operator displaces 
the phonon oscillation in the positive and negative di¬ 
rection. The parameter f n ^ determines the magnitude 
of the displacement for each mode. It is obvious that 
fnk = 9nk corresponds to the conventional polaron trans¬ 
formation, while f n k = 0 denotes no transformation. 

After the unitary transformation, the Hamiltonian of 
the system becomes 


Hs — ^ ^ (Cn T RnknCni (4) 

n 

where R n = \ n ,k{fn,k — ^gn,k) is the renormalization en¬ 
ergy which describes the shift of single-particle energy of 
electronic states induced by electron-phonon interaction. 
H a is same as H a since it commutes with the unitary op¬ 
erator U. While, the system-lead coupling Hamiltonian 
becomes 

H Sa = Y} Vk ° nC L CnXn + Hc -] ( 5 ) 

nk a 

after the transformation, i.e., the system-lead coupling is 
dressed by phonon displacement operator X n: where X n 
is defined as 


X n = exp 


^4 An,fc(4,fc b n ,k) 

k 


( 6 ) 


After the unitary transformation, the Hamiltonian of 
electron-phonon coupling becomes 


Hsb — ^ ' flre,fc4 t c n(4.fc 4” ^n,k)> (7) 

n,k 


where g n ,fc = 9n,k — fn,k is defined as the reduced 
electron-phonon coupling strength. g nt k indicates that 
the electron-phonon coupling is reduced by the unitary 
transformation. In summary, after the variational po¬ 
laron transformation, phonon’s influence on electrons is 
transferred into three terms: (1) the polaron shift energy 
R n , which represents the energy renormalization effect 
due to electron-phonon interaction; (2) in the dressed 
coupling between system and leads, the system-lead cou¬ 
pling is renormalized by the displacement operator X n ; 
(3) reduced electron-phonon coupling g n ,k, which is zero 
in the strong coupling limit and approaches g n k in the 
weak coupling limit. 

The variational transformation results in a new kind of 
system-lead coupling, which is mediated by the phonon 
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displacement operator. The displacement operator can 
be divided into two parts: 

X n = X n - (X n ) + (X n ) = B n + (X n }, 

i.e., the expectation value and corresponding fluctuation 
around its expectation. Assuming phonon is in thermal 
equilibrium, ( X n ) can be written as 


(Xn) = exp 


1 

2 


E A n,fe C0th (/3w„,fe/ 2 ) , 

k 


( 8 ) 


where /? is the inverse temperature. Consequently, Hs a 
can be separated into two parts, 

^5a=E^« C L^+H.C.], 

n,k a 

HL =Y,lVk a ncl a c n B n +K.c.]. (9) 

n,k a 

where V^ n = Vk an (X n ) is the renormalized system-lead 
coupling mediated by phonon. Thus, the total trans¬ 
formed Hamiltonian can also be grouped into two parts. 
One part is the mean-field (or non-interacting) part Ho 
and the other one is the interacting Hamiltonian Hi, 
where H 0 is defined as 


H 0 =H s + '52(H% a + H a ) + H B =He + Hs- (10) 

O' 

and Hi is 

Hi = H sb + E (11) 

a 

It is obvious that the polaron transformation (f n ,k = 
gn,k) makes the electron-phonon interaction term Hsb 
vanished. In contrast, the variational polaron transfor¬ 
mation employs f n k as parameters and the values of 
which are determined by the minimization of an upper 
bound of free energy. The minimization inherent to vari¬ 
ational approach allows to minimize the effect of Hi and 
obtain an optimal mean-field Hamiltonian Hq. Thus, 
given arbitrary electron-phonon coupling strength, the 
variational polaron transformation in principle can end 
up with small Hi. Consequently, 2PT can be employed to 
account the effect of Hi. Because of the single-particle 
nature of Ho and second-order perturbative treatment 
of Hi, variational polaron theory significantly reduces 
the computational complexity compared to the many- 
particle approaches. 


C. Minimization of the upper bound of free energy 

In this work, the Feynman-Bogoliubov upper bound of 
free energy 32 is chosen to determine the optimal value of 
f n k for the variational transformation. The upper bound 
of free energy is given as 


where (Hi)h 0 = Tr [Hie~° H °]. F u gives the upper bound 
of the true Free energy F, which is related to F by the 
inequality F < F u . It is obvious that the second term 
on the right hand right (RHS) of Eq.(12) is zero by the 
construction. Besides, because [Hb, H e ] = 0, F u can be 
further simplified and divided into two parts 

Fu=F B -±ln{tT[e-P S °]}, (13) 

where Fb is the free energy of phonon bath. Because Fg 
is not a function of f n k, it is neglected in the minimization 
procedure. There are two ways to minimize the upper 
bound of free energy with respect to f n k■ One way is 
to treat the system and leads as a whole system and the 
upper bound of free energy can be numerically calculated 
by the diagonalization of H e . The other approach is to 
employ the Linked cluster expansion method. 39 

The mean-field part of the transformation Hamilto¬ 
nian, Ho, is a function of the phonon-induced renor¬ 
malization parameters {R n , (X n )}. Consequently, the 
Feynman-Bogoliubov upper bound of free energy F u can 
be minimized with respect to {R n , (A' rl )}. Hence, the 
minimization condition can be written as 

dA s dA s dR n dA s d(X n ) 

dfn,k dRn dfn,k d(X n ) df n , k ' ^ ' 


Using the expressions for {R n , (X n )}, we can evaluate 
the derivatives of {R n , ( X n )} with respect to f n ,k analyt¬ 
ically. After simple algebra, the minimization condition 
gives the expression of the variational transformation pa¬ 
rameter /„ jfc as /„ jfc = g n ,kF n ,k, where 


Since {R n , ( X n )} are also functions of f n k, above equa¬ 
tion must be solved self-consistently. After the minimiza¬ 
tion procedure, the optimal choice of Hq is obtained, re¬ 
sulting in a small Hi. Hence, Hi can be treated by 2PT. 
Next, dissipative quantum transport formalism with vari¬ 
ational polaron transformation is presented. 


III. DISSIPATIVE QUANTUM TRANSPORT 
THEORY WITH VARIATIONAL POLARON 
TRANSFORMATION 


Non-equilibrium Green’s function (NEGF) formalism 
is employed to study the dissipative quantum transport. 
The key quantity in the NEGF method is the single¬ 
particle Green’s function, which is given by 

Gij(T, T ') =-*(T c c i (r)c](r , ))i/ 

=-i(T c Ci(T)Xi(T)c)(T')X](T'))H, (16) 


Fu = ~\n{tv[e-^]} + (Hi) So 


( 12 ) 


where H and H are the Hamiltonians before and after the 
variational polaron transformation, respectively, r and 
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t' are the time variables defined on the Keldysh contour 
and T c is the contour time-ordering operator. Eq.(16) 
determines the dynamics of the coupled electrons and 
phonons. We employ the following approximation to de¬ 
couple the electron and phonon dynamics 

Gij{T,T r ) = G ZJ (T,T')K ij {T,T'), (17) 

where 

Gij (t, t') =-i{T c Ci(r)c){T')) h 

^■(T,r')=(T c W(r)Xj(r')) S . (18) 

The decoupling in Eq.(17) is inherent in the Born- 
Oppenheimer approximation. Even the decoupling 
approximation is made, there is still correlation be¬ 
tween electron and phonon if self-consistent procedure 
is operated, 26 which is similar to the diagram dressing 
process in the standard many-body perturbation theory. 

With the decoupling approximation, the key quantity 
becomes G(r, r') which is associated with H. As illus¬ 
trated in Sec. II, the variational polaron transformed 
Hamiltonian can be divided into two parts: Hq and 
Hj. Hq is mean-field part of the Hamiltonian, the quan¬ 
tum transport problem associated with Hq can be easily 
solved as demonstrated by previous publications. 40 While 
the correction of Hi to the quantum transport can be 
taken into consideration by the perturbation theory since 
the variational polaron transformation in principle ends 
up with small Hi. Firstly, zero-order Green’s function 
Gq(t,t') is obtained by using H 0 , 


where E/(t, t') is the self-energy to Green’s function at 
the second-order with respect to Hi and Go(r, r') is the 
Green’s function corresponding to Hq. Since Hi contains 
two parts, H : Sa and Hsb , E/(t, t') can also be divided 
into two parts accordingly, 

E/(t, t’) = E ep (r, t') + E E/«(t, t'). 

a 

Within the 2PT, E i(t,t') is derived as 

S/cw (t, t') =E V£ ai V kaj g ka (r, r') 

k a 

x {T C Bj {t')B\(t))ii o . (22) 

The definition of E/ a is very similar to Eq.(20) except 
the phonon renormalization factor (X,) is replaced by the 
correlation function of B^. According to the definition of 
Bi , the self-energy E Q / can be rewritten as 

Sj a>i j(r, t') =T, a i j (r, r')[Kji - (X i ){X j )] 

=S a ,y(r, T')Kji(T',T) - E 0 a ,ij(T, r')(23) 

where EoaG, r') is the phonon dressed self-energy, which 
is used to evaluate the zero-order Green’s function 
Go(r, t') as described previously. 

The derivation of E ep (r, r) is same as that in Ref. 6. 
The only difference is that the electron-phonon coupling 
constant (g n , k ) is replaced by the reduced one ( g n ,k )• The 
expression of E ep (r, r') is 


G 0 (t,t') =go{r,r') + J dn J dT 2 go(T,Ti) 


E 0q (ti, r 2 ) 


Go(T2,r'), (19) 


where go(i~, t') is the electron Green’s function in absence 
of system-lead coupling. The effect of system-lead cou¬ 
pling is represented by the self-energy 


Soa,«(r, r’) =Y, V k * ai V kaj g ka (r, r')(X l )(X J ) 
k a 

=E a , lJ (r,r')(X i )(X i ), (20) 

where g ka {T,T') is the surface Green’s function of lead 
a and E a (r, r') is the lead self-energy in absence of 
phonon renormalization. Next, G(r, r') is derived at 
second-order expansion with respect to U/, starting from 
Gq(t, t'). 


A. Second-order perturbation theory 


Within the 2PT theory, the Green’s function can be 
expressed as 

G(t,t') = G 0 (t,t')+ J dr i J dT 2 G 0 (T,Ti) 

xE/(ti,t 2 )Go(t 2 ,t'), (21) 


S ep (r, t') = ^2 9 n,kD Q {T , r')G 0 (r, T r )g n}k , (24) 

n,k 

where Do(r, r') is the bare phonon Green’s function. 

In practical implementation, the Green’s function and 
self-energies defined on the Keldysh contour are projected 
on real axis, resulting in equations for retarded/advanced 
and lesser/greater Green’s functions and self-energies. In 
this work, we focuses on the steady state. In steady 
state, all quantities (two-time functions) are only depen¬ 
dent on the time difference instead of two times. Hence, 
Fourier transformation is employed to transfer the quan¬ 
tities from real-time domain to energy space. Since Hq 
is the mean-field part of the total Hamiltonian, the cor¬ 
responding Green’s functions can be easily evaluated, 

G r 0 (E) =[E + ir] -HQ~J2 Z r o a (E)}~\ 

a 

G<(E) =G5(E)[^E< q (E)]G 0 “(E), (25) 

O' 

where rj —> 0 + . E(j a and Eg a are phonon dressed retarded 
and lesser self-energies of lead a, respectively, which are 
written as 

=E v*Vk aj g r k / a < (E)(x i )(x j ) 

ka. 

=E^(£)(X i >(X i ) 


(26) 
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Defining the unperturbated line-width function as 
T a (E) = — 23[E^(.E)], where 3 denotes the imaginary 
part. Thus, the phonon dressed line-width function be¬ 
comes f a ,ij(E) = T ai j j(E)(Xi)(Xj). The dressing fac¬ 
tor indicates the localization of electron induced by the 
electron-phonon interaction. In term of line-width func¬ 
tion, the corresponding lesser self-energy can be rewritten 
as Y,<(E) = if a {E)T a {E) according to the fluctuation- 
dissipation law, where f a (E) is the Fermi-Dirac distribu¬ 
tion function of lead a with chemical potential \x a . 

After obtaining the zeroth-order Green’s functions, 
G r J < {E), are evaluated, the Green’s function with Hi 
is obtained from the 2PT theory as 

G r (E) =Gm+G r 0 (E)^(E)G r o(E), 

G<{E) =G<(E) + G r 0 (E)X<(E)G<(E). (27) 

where E^' /< (£') is the self-energy to Green’s function at 
the second-order with respect to Hi. As shown previ¬ 
ously, E’/^.E) contains two parts, 


Retarded and lesser phonon self-energies within the 2PT 
are 6 - 41 ' 42 


(E) =]T g q N±G< {E ± Uq )g q , 

9,± 

K P { E ) =E 9 q[N*G r 0 {E ± w,) ± G<(E t Uq )]g q { 28) 

9,± 


where N± = N q + ^ ± i. N q is the phonon occupation 
number of mode g, which is determined by Bose-Einstein 
distribution function. 

The derivation of requires the knowledge of 

the displacement correlation function AT(r, t'), which is 
nontrivial. K(t,t') in principle depends on the phonon 
Green’s function which is coupled with electronic Green’s 
function via its self-energy. 26-28 Therefore, self-consistent 
calculation of phonon and electron Green’s function is 
required. However numerical implementation of self- 
consistent calculation is non-trivial and very computa¬ 
tional demanding. In practise, we assume the phonon is 
in equilibrium and undressed by the electron in this work. 
The influence of electron to the phonon can be intro¬ 
duced through a phenomenological rate equation includ¬ 
ing the renormalization, damping and heating effect. 43-45 
With the assumption that phonon is in the equilibrium 
and undressed by electron, the displacement correlation 
function can be rewritten in a simple form. 24 ’ 26,39 The 
lesser projection of displacement correlation function is 
expressed as 


= n (e-^ , -.«)exp{A if A jt x 

9=1 l 

TV,e+ (TV, + i) e < "«( t -*')] } 


where N q is the occupation number for the gth phonon 
mode. K<j(t,t') can be decomposed as 


M 




9=1 


EA' 


,ip q ui q (t-t') 


Pi 


= E rry:y.. L ~vi> T »<«-<'> 






(30) 


where both p and cu are row vectors, p T u> = 

And Lfj = Lfj ■ ■ ■ L ™, where L 1 ’’’ is modified Bessel 
function 


L P? =e ~ Xiq+ 2 Xjq ( 2N q +l) e p q u q p/2 


xl. 


Pi 


2\ iq \j q ^ TV, (TV, + 1) 


(31) 


and I Pq is the p,th order Bessel function. 

The greater projection of displacement correlation 
function is = (Xi(t)Xt(f)) = f')] t . It 

can be verified that K < (t,t') ~ K > (t,t') in the high- 
temperature limit since TV, ~ TV, + 1 in this limit. Actu¬ 
ally, many works even adopt K < {t 1 t’) ~ as an 

approximation in low temperature, which is regarded as 
neglecting the Fermi Sea. 22-25,46 

With K < ’ > (t,t') known, the lesser (greater) self¬ 
energy in energy space can be obtained from the Fourier 
transformation of E (t, t'), 

^ 3 {e) L pKiA E - p T “) - s; Uij(E), 

p 

^j(E) =E UjiKiA E + P T ^ - ^a,ij(E). (32) 

p 


B. Observable of interest 

The steady-state current through the junction can be 
expressed in Meir-Wingreen formula 40 

= T / f T «< £ > < 33 > 

where the factor 2 accounts for the spin degeneracy and 
T a is the transmission coefficient through lead a , which is 
expressed in terms of Green’s functions and self-energies 
within NEGF formalism 

T a {E) = Tr [E <(E)G>(E) - E >(E)G<(E)] . (34) 

As shown in Eq.(17), lesser/greater projection of Green’s 
functions are 


m 
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Utilizing the expansion of K?j(t,t'), Fourier transforma¬ 
tion of above equations gives the Green’s functions as 

G>(E) =Y J Ll l Gf ] (E-p T u;), 

P 

G<(E) =£ L*G<(E + p T u>). (35) 

p 

Hence, transmission coefficient and current can be ob¬ 
tained through substituting Eq.(35) into Eq.(34). 

In short, the computational procedures are summa¬ 
rized as follows: (1) Firstly, bare Green’s functions 
Gq^(E) in variational polaron picture are calculated; (2) 
Secondly, Green’s functions G r ’>(E) at the second-order 
with respect to Hi are evaluated; (3) Thirdly, Green’s 
functions in the H picture are obtained through Eqs.(35); 
(4) Finally, physical quantities, such as density of states 
(DOS) and currents, are obtained from the Green’s func¬ 
tions and self-energies. 


value / of variational polaron transformation is searched 
for different phonon frequency and electron-phonon cou¬ 
pling strength. The phonon frequency is chosen in three 
different regimes: uj > t, uj = t and uj < t. The op¬ 
timal value / of variational polaron transformation for 
different electron-phonon coupling constant g is shown in 
Fig. 2. As it is expected, different electron-phonon cou¬ 
pling strength corresponds to different optimal value of 
/ in the variational transformation. Besides, / increases 
monotonously with electron-phonon coupling strength g. 
When g —> 0, f / g ratio approaches a small value which 
is affected by the phonon energy. In the g —0 limit, 
//g ratio is small, which is close to no polaron transfor¬ 
mation. While in the strong electron-phonon coupling 
limit ( g —> oo), f/g ratio approaches to 1, which corre¬ 
sponds to full polaron transformation. In the intermedi¬ 
ate regime, variational polaron transformation ends up 
with f/g < 1, which denotes the partial polaron trans¬ 
formation. 


IV. NUMERICAL EXAMPLES 

The variational polaron theory for dissipative quantum 
transport is employed to study the transport through on 
a monoatomic chain. The system consists of a single site. 
The two leads are modeled by semi-infinite linear chains 
as illustrated in Fig. 1. The onsite energy of the system is 
denoted by eg. The onsite energy and hopping element of 
leads are e and f, respectively. The system-lead coupling 
constant is denoted by V a . In this work, it is assumed 
that the system is coupled to the two leads equally, i.e., 
Vl = Vr = v. The line-width function of the semi¬ 
infinite tight-binding lead can be evaluated analytically 4 ' 



g/© 


r«(£) = 


V ^ 2 


(E-p a y 

t 2 


where p a is the chemical potential of lead a. 


FIG. 2. Optimal value / of variational polaron transforma¬ 
tion against different electron-phonon coupling constant g. 
Parameters are: The system-lead coupling and hopping in 
leads are v = t = 0.1 eV; T = 300 K. 


Left lead System Right lead 


t 

V V 

t 

• • • • 

~1T 

• • • • 

c 

£s 

£ 


FIG. 1. Illustration of a single site model coupled to two 
leads. It is assumed that the system coupled to the two lead 
equally, which is denoted by v. The onsite energy and hopping 
element of leads are e and t, respectively. 

In this work, the system-lead coupling is set to be equal 
to the hopping in the leads, i.e., v = t. The onsite energy 
is set to be 0 eV and hopping is set as t = 0.1 eV. And the 
Fermi energy of the leads is 0 eV, which corresponds to 
half-filled on each site. Temperature is set as T = 300 K. 
In this work, only one phonon mode is considered to be 
coupled to the system for simplification. The optimal 


It is also noticed from Fig. 2 that there is a shape tran¬ 
sition of / in a certain regime of g if uj < t. Larger 
results in the steeper transition. Finally, / shows discon¬ 
tinuity against —. This discontinuity appears because 
different local minimum becomes the global minimum at 
certain parameter regime. Taking uj = 0.05 eV as ex¬ 
ample, the free energies against f/g ratio with electron- 
phonon coupling strength around 0.1 eV are plotted in 
Fig. 3. As it shown by Fig. 3, two local minimums start to 
emerge at g ~ 0.1 eV. The local minimums of smaller and 
larger f/g become the global minimum for g = 0.098 eV 
and g = 0.102 eV, respectively. The discontinuity of / 
may cause unphysical phase transition from delocalized 
to localized regime. 33,48,49 For multiple sites systems, the 
free energy can be more complex and more local mini¬ 
mums may emerge, thus there is larger possibility of dis¬ 
continuity. However, even the discontinuity appears, the 
variational transformation is still likely to work better 
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than the untransformed scheme or polaron transforma¬ 
tion. Moreover, the discontinuity of / may be removed 
by employing certain ansatz. For instance, a ground state 
ansatz has been proposed recently to get rid of the dis¬ 
continuity in the variational polaron. 49,50 Similar ansatz 
may be proposed to remove the discontinuity problem 
in the quantum transport system, which is the scope of 
further improvement of the method. 



FIG. 3. Free energy as function of variational transformation 
parameter / for different electron-phonon coupling strength 
g. The phonon energy is 0.05 eV. When g is around 0.10 ~ 
0.102 eV, different local minimum may become the global 
minimum, resulting in abrupt transition of /. 

After the variational polaron transformation is found, 
the quantum transport properties of the system, such 
as DOS and current-voltage characteristics, at different 
electron-phonon coupling regime can be studied by the 
method developed in Sec. III. Here, three sets of param¬ 
eters are studied, ranging from weak to strong electron- 
phonon coupling regime. The parameters as well as opti¬ 
mal factor of variational polaron transformation are sum¬ 
marized in Table I. Model A and C are in the weak 
and strong coupling regime, respectively, and Model B 
is in the moderate coupling regime. The optimal values 
of variational polaron transformation are 0.17, 0.54 and 
0.90. Here, the g/u ratio is set as a constant for the three 
models. As shown by Fig. 2, large g/u ratio (g/u 1) 
makes the system in polaron regime while small g/u ratio 
(g/u « 1) makes the system in weak or moderate cou¬ 
pling regime. Hence, g/u = 0.8 is chosen to make sure 
that the system varies from the weak to strong coupling 
by tuning the electron-phonon coupling constant. 


TABLE I. Parameters of the Models. 


Model 

> 

CD 

3 

9 (eV) 

f/g 

A 

0.05 

0.04 

0.17 

B 

0.15 

0.12 

0.54 

C 

0.40 

0.32 

0.90 



E(eV) 

FIG. 4. DOS of the system with different electron-phonon 
coupling constant: (a) u = 0.05 eV and g = 0.04 eV; The 
electron-phonon interaction is in the weak coupling regime. 
DOS of variational polaron is same as that of 2PT, while 
polaron theory overestimates the phonon side-bands, (b) 
u = 0.15 eV and g = 0.12 eV; In this moderate electron- 
phonon coupling regime, 2PT starts to fail and polaron still 
overestimates the phonon side-bands, (c) to = 0.40 eV and 
g = 0.32 eV. The system is in the strong coupling regime. 
Variational polaron theory gives similar results as polaron the¬ 
ory does. While 2PT is not applicable to this regime, it even 
gives unphysical DOS. See text for other parameters. 

In presence of electron-phonon coupling, electronic 
structure of the system is affected at different levels de¬ 
pending on the electron-phonon coupling strength. Con¬ 
sequently, DOS of the system is affected, such effects in¬ 
clude energy shift of peak, broadening of the spectral 
and emergence of phonon side-bands. Fig. 4 plots the 
DOS of the system with different electron-phonon cou¬ 
pling strength. 

(a) In model A, the electron-phonon coupling strength 
















is small. Phonon has limited effects on the electronic 
structure of the system. The phonon-induced energy 
shift is about 50 meV. The DOS of variational polaron 
transformation is close to that of 2PT. Moreover, the 
DOS of both variational polaron theory and 2PT is al¬ 
most the same as that of non-interacting system except 
the phonon-induced shift of the peak. No obvious phonon 
side-bands are observed in the DOS of variational po¬ 
laron theory and 2PT. However, polaron theory over¬ 
estimates the energy shift and phonon side-bands. Be¬ 
cause polaron theory overestimates the phonon dressed 
system-lead coupling, resulting in stronger localization 
effect. Consequently, main peak of DOS is sharped and 
phonon side-band appears as shown in Fig. 4(a). 

(b) In the moderate coupling regime, phonon has 
larger effect on the electronic structure of the system and 
phonon side-bands start to appear due to the increased 
phonon-induced electron localization. In this coupling 
regime, the DOS of variational polaron theory starts to 
deviate from that of 2PT as shown by Fig. 4(b). Due to 
the principle of 2PT, only two phonon side-bands can be 
observed in the DOS. However, higher order corrections 
become more and more significant in the moderate cou¬ 
pling regime. As a result, 2PT starts to fail to describe 
the renormalization effect induced by phonon. Besides 
the phonon side-bands, the energy shift due to electron- 
phonon interaction is underestimated slightly by 2PT. In 
contrast, higher order phonon side-bands can be evalu¬ 
ated by both polaron theory and variational polaron the¬ 
ory. But, polaron theory still overestimates the phonon 
dressed system-lead coupling compared with variational 
polaron theory. Consequently, the phonon side-bands are 
much more significant in the DOS. While, variational 
polaron theory neither underestimates nor overestimates 
the renormalization effect. The DOS of variational po¬ 
laron lies between those of 2PT and polaron. 

(c) In the strong coupling limit, phonon induces pro¬ 
nounced side-bands effect and polaron is formed. Con¬ 
sequently, the variational polaron transformation ap¬ 
proaches the polaron theory. As shown by Fig. 4(c), 
DOS of variational polaron theory is almost the same 
as that of polaron theory since the optimal value of 
f/g is 0.9, which is close to 1, i.e., variational polaron 
theory approaches the polaron theory. In the contrast, 
2PT greatly underestimates the phonon side-bands and 
overestimates the phonon-induced energy shift. The 
pronounced phonon side-bands indicate that electron is 
localized by the phonon. The phonon side-bands is 
also responsible for the phonon blockade in the current- 
voltage characteristics, which is represented by steps in 
the current-voltage characteristics or peaks in the differ¬ 
ential conductance. 23-26 


V. SUMMARY 

In this work, a dissipative quantum transport the¬ 
ory with electron-phonon coupling at arbitrary param¬ 


eter regime is developed by employing the variational 
polaron transformation. Variational polaron transforma¬ 
tion provides an optimal set of parameters for the uni¬ 
tary transformation, which is determined by the min¬ 
imization of the Feynman-Bogoliubov upper bound of 
free energy. After the minimization, variational polaron 
transformation is able to end up with an optimal mean- 
field Hamiltonian and a small perturbation. The mean- 
field part Hamiltonian, Hq, contains the phonon-induced 
energy renormalization and phonon dressed system-lead 
coupling. Thanks to the variational polaron transfor¬ 
mation, the many-body effect of electron-phonon inter¬ 
action is transformed into the small interacting part of 
the Hamiltonian ( Hj ), which validates the 2PT treat¬ 
ment. Upon the variational transformation, a quantum 
transport theory with electron-phonon interaction at ar¬ 
bitrary parameter regime is established within NEGF for¬ 
malism. Following the quantum transport theory of non¬ 
interaction system, the mean-field part is treated exactly. 
Hi is taken into account through the self-energies which 
are in second-order of Hi 

Numerical examples on a monoatomic chain demon¬ 
strate the validity of the variational transformation. 
The optimal value of variational transformation increases 
monotonously with electron-phonon coupling strength, 
which connects the 2PT in weak coupling regime and 
polaron theory in strong coupling limit. Comparison 
of DOS in different parameter regimes indicates the ap¬ 
plicability of variational polaron transformation: (a) In 
weak coupling limit, polaron theory overestimates the 
phonon-induced energy shift and side-bands; (b) In the 
strong coupling limit, 2PT underestimates the renormal¬ 
ization effect induced by phonon and even gives unphys¬ 
ical DOS; (c) Variational polaron transformation natu¬ 
rally connects the 2PT in the weak coupling limit and 
polaron theory in the strong coupling limit by making use 
of the optimal transformation. More numerical bench¬ 
mark comparison with other numerically exact methods, 
such as QUAPI and MCTDH, will be considered in the 
future works. Moreover, the sudden change of the opti¬ 
mal value of variational transformation may be solved by 
employing certain ansatz, 49,50 which is also the scope of 
future improvement of the method. 
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